function   [Winv,WFactor,WLogDet]=svdInverse(W,tolZero)
% function [Winv,WFactor,WLogDet]=svdInverse(W,tolZero)
% 
% For a symmetric matrix W Compute the inverse and, if requested, 
% a Cholesky type factorization and the log determinant using the SVD. 
% Small singular values are ommitted. 
% 
%% *Inputs* 
%  W            [n,n] symmetric matrix 
%  tolzero      (optional) maximum tolerance for small singular values 
%                default is 1e-12 
% 
%% *Outputs* 
% Winv          inverse of W 
% WFactor       (if two outputs requested) cholesky-type factor such that W=WFactor'*WFactor 
%               WFactor is NOT upper triangular as if one used chol(W)
% WLogDet       (if three outputs requested) log( det(W) ) 
% 
% Alejandro Justiniano 
%==========================================================================
if nargin < 2 || isempty(tolZero)==true 
    tolZero=1e-12; 
end 
W=0.5*(W+W'); 

%% 1  SVD(W)
% PP*DD*PPinv'=W  Notice the transpose 
% PPinv=inv(PP)' 
% PPinv'=inv(PP); 
% PP=PPinv
[PP,DD,PPinv]=svd(W); 
%comparemat( PP , PPinv ); 
%comparemat( PP\eye(length(PP)), PPinv' ); 


%% 2. Truncate small singular values 
firstZero=min(find(diag(DD) < tolZero));
if isempty(firstZero)==false 
    warning('Truncating singular values for inversion'); 
    PP=PP(:,1:firstZero-1); 
    PPinv=PPinv(:,1:firstZero-1); 
    DD=DD(1:firstZero-1,1:firstZero-1); 
end 

%% 3. Inverse 
Winv=PP*diag( 1./diag(DD) )*PPinv'; 
%comparemat(Winv,W\(eye(length(DD)))); 
if nargout < 2; return; end 

%% 4. Cholesky type factorization
%     NOTE: the matrix is NOT upper triangular 
WFactor=diag(sqrt(diag(DD)))*PPinv'; 
%comparemat(WFactor'*WFactor,W); 

if nargout < 3; return; end 
%% 5. determinant 
WLogDet=sum(log(diag(DD))); 
%comparemat(WLogDet,log(det(W))); 